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I .  Introduction 

The  purpose  of  this  volume  will  be  to  provide  an 
evaluation  of  the  various  algorithmic  options  inherent  to 
ADINA  (1977) .  The  main  thrust  of  the  work  will  be  to  esta¬ 
blish  the  inherent  characteristics  of  the  static,  dynamic 
and  eigenvalue  extraction  solution  branches.  This  will 
be  considered  in  two  steps: 

i)  Overview  framework  of  algorithmic  options, 
and 

ii)  Establish  computational  capabilities  and  short¬ 
comings.  This  is  achieved  through  extensive 
benchmarking. 

The  benchmarking  utilized  will  itself  serve  two  basic 
purposes  namely: 

1)  Establish  the  algorithmic  sensitivities,  con¬ 
vergence  characteristics  and  typical  solution 
failures  generic  to  ADINA  and  secondarily; 

2)  Perform  code  check  out. 

For  the  static  and  dynamic  branches,  benchmarks  will 
be  used  to  quantify  the  convergence,  stability  and  artificial 
propagation  characteristics.  In  this  direction,  the  results 
of  several  numerical  experiments  will  be  employed  to  check 
out  the  handling  of: 

a)  Kinematic  nonlinearities; 

b)  Material  nonlinearities  and; 

c)  Combined  geometric  and  material  nonlinearities. 


Special  emphasis  will  be  given  to  ascertaining  the  operating 
characteristics  in  the  presence  of  softening,  hardening, 
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elastic/plastic  and  load/unload  situations.  In  this  con¬ 
text,  the  major  purpose  of  the  benchmarking  will  be  to 
establish  the  sensitivity/pathological  behavior  of  the 
various  generic  nonlinear  algorithms  to  changes  in: 

a)  Convergence  criteria 

b)  Time  step  size 

c)  Load  increment  size  and 

d)  Material  models  etc. 

e)  Order  of  integration  of  element  stiffnesses 

To  enable  the  evaluation  of  the  various  inherent 
anomalous  behaviorial  pathologies,  the  strain  energy  will 
be  monitored  during  the  various  stages  of  calculation.  This 
includes  both  local  and  global  evaluations.  Such  an  approach 
will  allow  for  the  monitoring  of  the  initiation  characteristics 
of  pathological  behavior  which  is  either  local  or  global  in 
nature . 

In  addition  to  testing  the  purely  transient  phase,  the 
eigenvalue/vector  extraction  routine  will  also  be  evaluated. 

The  benchmarks  used  in  this  phase  of  evaluation  will  establish 
such  factors  as  convergence  properties,  eigenvalue/ vector 
deterioration,  multiplicity  and  separability  etc. 

As  a  natural  outgrowth  of  the  foregoing  approach, 

various  aspects  of  the  configuration  control  potential  inhe- 

% 

rent  to  ADINA  will  be  discussed.  This  will  include  such 
features  as: 

i)  Equation  solver  options  (eigenvalue  extraction 
routines,  integration  schemes  as  well  as  various 
nonlinear  algorithms)  and; 
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ii)  Equation  solver  parameters  (time  steps,  load 
increments,  tolerance  limits  iteration  counts 
etc. ) 

Since  many  of  these  features  are  of  significant  importance 
to  nonlinear  codes,  each  of  the  foregoing  configuration 
control  options  will  be  examined  in  detail. 

Based  on  the  foregoing,  the  overall  report  will  be 


organized 

in  four  main  sections  including: 

i) 

An  overview  of  the  basic  framework  of 

linear  algorithmic  solution  options; 

the  non- 

ii) 

A  detailed  benchmarked  discussion  of 

capabilities  and  shortcomings; 

algorithmic 

iii) 

A  discussion  of  suggested  algorithmic 

and  lastly; 

improvements 

iv) 

A  summary  of  results. 

Appendices  which  include  equations  for  calculating  the  global 
and  local  strain  energies,  code  checkout  benchmarks  and  program 
listings  of  code  modifications  employed  in  the  study  and  lastly 
sample  problem  data  input  echoes. 
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II.  Overview  of  Framework  of  Algorithmic  Solution  Options 

As  noted  earlier,  this  section  will  consider  the  numerical 
framework  of  the  various  algorithmic  options  available  to  ADINA. 
In  particular,  this  will  include  discussions  of  the  generic/ 
static  nonlinear  equation  solver,  the  dynamic  and 
eigenvalue  solution  branches  /  the  convergence  schemes,  user 
I/O  options,  adaptive  strategies  and  programmed  error  stops. 

II.  1  Generic/Static  Nonlinear  Equation  Solver 

Based  on  the  theoretical  framework  of  the  governing  field 
equations  outlined  in  Part  I  it  follows  that  the  Modified 

Newton  Raphson  (MNR)  algorithm  is  the  generic  ADINA  nonlinear 
equation  solver.  In  particular,  it  is  employed  both  for  the 
static  and  dynamic  branches  of  the  code.  The  algorithm  is 
encoded  with  three  main  options  namely: 

i)  Load  incrementation  with  iteration  and  reformation 
(only  once) ; 

ii)  Load  incrementation  with  iteration  but  no  reformation 
and  lastly; 

iii)  Load  incrementation  and  reformation. 

The  overall  logic  flow  associated  with  such  options  is  defined 
in  Fig.  1.  In  particular  the  reformation  loop  in  ADINA  has 
been  specialized  so  as  to  provide  for  nonlinear  element  stiff¬ 
ness  updating  only  at  the  beginning  of  a  new  load  or  time  step. 
From  that  point  on,  iteration  proceeds  without  reformation. 
Figures  (2,  3,  and  4)  present  a  graphical  representation  of 
the  flow  of  calculation  associated  with  such  options.  Note 
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these  options  have  been  programmed  for  both  the  static  and 
dynamic  branches  of  ADINA. 

The  convergence  checks  associated  with  the  MNR  algorithm 
are  two- fold,  namely, 

1)  A  norm  test  of  successive  nodal  increment  vectors 
as  represented  by 

I i £i+l ^  I  —  tolerance 

2)  A  norm  test  of  successive  out  of  balance  loads  as 
represented  by 

I lR±l l  >  I  I Ri+i i I 

The  first  test  determines  whether  the  solution  has  converged. 
The  second  assures  that  the  tendency  of  convergence  is 
maintained. 

Note  while  the  iteration  counter  limit  can  be  superceded 
by  a  user  I/O  option,  the  out  of  balance  and  negative  pivot 
stops  do  not  allow  for  direct  user  intervention.  This  appears 
appropriate  since  such  events  are  usually  a  natural  outgrowth 
of  some  form  of  solution  divergence.  In  ADINA,  typically -such 
pathologies  occur  under  the  following  broad  category  of 
situations  namely: 

1)  In  the  neighborhood  of  buckling  points  (Fig.  II.5)  wherein 
there  is  either  a  significant  change  or  a  complete  reversal 
of  structural  stiffness; 

2)  In  the  neighborhood  of  plastic  unloading  zones.  Fig.  n.6; 


1 1  -  6  Typical  Load  Deflection  Curve  for  E 1  a s t i c- P 1  a s t i c  Load-Unload 
Problem 


3)  Use  of  excessive  load  step  size  and; 

4)  Use  of  excessive  time  step  size. 

Typically,  in  the  vicinity  of  buckling  points  either 
the  out  of  balance  load  or  negative  pivot  stop  are  encountered. 
As  seen  from  the  scenarios  depicted  in  Figs.  11.7  and  r,  such 
stops  are  a  direct  outgrowth  of  the  use  of  the  MNP.  operator 
which  is  driven  by  the  "tangent  hvrer-plane"  of  the  solution 
surface.  Note  the  tangent  plane  in  the  neighborhood  of  buck¬ 
ling  points  can  become  quite  shallow  with  possible  "slope" 
reversals.  Because  of  this,  as  currently  programmed,  there  is 
no  possibility  to  analyze  the  post  buckling  range  through  the 
user  I/O  intervention  capabilities  currently  available  in  ADINA. 
What  is  obviously  necessary  is  the  incorporation  of  heuristic 
(adaptive)  programming  to  circumvent  the  difficulties  in  such 
problems.  This  will  be  discussed  in  more  detail  later  in 
Section  IV. 

The  reason  for  solution  blow  ups  in  the  neighborhood  of 
plastic  unloading  zones  is  essentially  a  direct  outgrowth  of 
the  modified  NR  algorithm  employed  by  ADINA.  As  reformation 
occurs  only  at  the  beginning  of  a  load  step,  if  some  unloading 
occurs  during  the  iteration  phase,  then  it  is  not  properly 
accounted  for.  As  will  be  seen  later,  this  obviously  leads  to 
solution  drift  and  typically  out  of  balance  loads. 

While  the  use  of  excessive  load  step  sizes  can  lead  to 
solution  stops  generated  by  any  of  the  checks,  typically  the 
out  of  balance  stop  is  most  often  encountered.  Generally, 


Buckling  Zone 
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excessively  large  steps  tend  to  lead  to  solution  divergence. 
Usually  this  is  due  to  the  fact  that  large  load  steps  may  excite 
significant  changes  in  the  stiffness  of  a  given  structure.  If 
these  changes  occur  within  a  given  load  step,  the  modified  INR 
employed  by  ADINA  will  miss  them.  Hence  since  the  hyper  surface 
representing  the  solution  is  extremely  complex,  the  algorithm 
will  usually  drift  aimlessly  over  the  surface  until  typically 
out  of  balance  loads  are  excited.  In  a  similar  manner,  excess¬ 
ive  time  step  sizes  may  also  lead  to  solution  divergence  caused 
by  the  large  changes  occurring  during  a  given  interval  of  time. 

As  was  noted  earlier,  the  1977  version  of  ADINA  does  not 
possess  any  heuristic  programming  which  can  automatically  cir¬ 
cumvent  convergence  difficulties.  The  only  recourse  open  to  the 
ADINA  user  is  to  employ: 

1)  The  restart  option  and/or; 

2)  The  time  block  option. 

With  these  options,  it  is  possible  to  either  manually  reset  the 
algorithmic  option  and  restart  or  preset  the  algorithmic  choice 
at  a  predetermined  point  in  time.  Generally  the  scenario  used 
follows  the  sequence  of  events  defined  below: 
i)  Restart  option 

i.  Store  all  field  quantitites  generated  on  tape; 

ii.  When  error  stop  is  encountered,  manually  restart 
with  new  algorithmic  choice  and; 


iii.  Continue  process  to  successful  completion. 
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2)  The  Block  Option 

i.  Having  employed  the  restart  option  to  obtain  a 
successful  solution,  points  of  unloading  can  ob¬ 
viously  be  aetermined; 

ii.  With  such  information  available,  the  appropriate 
algorithm  can  be  selected  for  an  estimated  time 
block. 

iii.  Once  all  the  time  blocks  are  created,  ADINA  will 
automatically  use  the  preselected  algorithm  for  a 
given  time  interval. 

While  the  foregoing  options  are  certainly  an  improvement 
over  the  straight  primitive  algorithmic  choices,  it  presupposes 
user  knowledge  of  the  various  solution  sensitivities  which  for 
nonlinear  problems  are  both  structure  and  loading  dependent. 
That  is,  the  solution  sensitivities  of  a  given  nonlinear  struc¬ 
ture  may  shift  from  load  to  load.  Furthermore,  for  "one  of 
the  kind"  problems  for  which  nonlinear  codes  are  oftentimes 
used,  such  knowledge  can  be  obtained  only  after  sufficient 
parametric  studies  have  been  completed.  As  noted  earlier,  such 
an  approach  can  oftentimes  require  significant  manpower  allo¬ 
cations  as  well  as  extensive  amounts  of  computer  time. 

II. 2  Dynamic  Equation  Solver 

As  noted  in  Part  I,  for  time  dependent  problems  numerical 


integration  is  required  to  obtain  a  solution.  In  ADINA,  this 
is  achieved  via  the  use  of  either  implicit  or  explicit  opera¬ 
tors  namely: 
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1)  Implicit  operators: 

i)  Wilson  0  method^ 

■  •  r-i 

u)  Newmark  method  ^ 

2)  Explicit:  Central  Differences^ 

To  trace  through  the  algorithmic  flow  of  the  transient  branch, 

the  use  of  the  Wilson  operator  will  be  considered.  Its  form 
,  [2  ] 

is  given  by: 


t+At  .. 

1  ^ 

.. 

t+0 At  .. 

AY  = 

'X, 

u-  |) 

AY 

+  e 

AY 

(ii. l) 

t+At 

AY  = 

t  • 

AY  + 

At 

( AY 

t+At  .. 

+  AY) 

(II. 2) 

'X, 

'X/ 

2 

V 

-x. 

t+At 

AY  * 

'X, 

fcAY  + 

O/ 

<t 

rt 

AY  + 

V 

,  ...  2  t+At  ..  .  .. 

(^r  (  AY  +  2fcAY) 

(II. 3) 

Employing  the  operators  defined  by  (II. 1  -  3),  the 
following  sequence  of  operations  is  necessary  to  develop  the 
solution  for  one  time  step  namely: 

1.  At  beginning  of  new  time  step  update  structural  stiffness; 

2.  Form  effective  load; 

3.  Solve  for  incremental  displacements; 

4.  If  required,  iterate  dynamic  algorithm  and; 

5.  Upon  convergence,  establish  final  fields  for  beginning  of 
next  time  step  (displacements,  velocities,  accelerations, 
stresses,  strains,  etc.). 

Since  the  generic  nonlinear  equation  solver  employs  a  modified 
form  of  NR  algorithm,  the  iteration  process  designated  by  Step  4 


involves  no  reformation  of  the  stiffness.  As  with  the  static 
solution  branch,  the  convergence  check  associated  with  the 
iteration  process  involves  the  use  of  norm  tests  on  the  dis¬ 
placements  and  out  of  balance  loads.  In  short,  apart  from 
the  capability  of  choosing  either  implicit  or  explicit  opera¬ 
tors,  all  the  iteration  options,  convergence  checks,  error 
stops,  user  intervention  options  (restart,  time  blocks)  etc. 
are  the  same  as  those  of  the  static  branch. 

II.  3  Eigenvalue/vector  Algorithm 

Since  the  calculation  of  the  frequencies  and  associated 
mode  shapes  is  important  in  modal  analysis,  many  general  pur¬ 
pose  codes  have  been  provided  with  an  extraction  algorithm.  In 
nonlinear  codes  such  as  ADINA,  the  main  importance  of  the  eigen 
value/vector  algorithm  lies  in  the  fact  that  it  can  be  used  to 
establish  the  requisite  time  increment  for  dynamic  analysis 
problems . 

In  ADINA  the  algorithm  considers  the  generalized  eigen¬ 
value  problem  defined  by 

[K]  Y  =  co2  [M] V  (II. 4) 

where 

[K]  -  is  the  tangent  stiffness  matrix  at  time  zero 

[M]  -  is  the  mass  matrix 

V  -  the  modal  vector 
% 

co  -  the  system  frequency 

Because  of  the  generality  of  programming,  the  mass  matrices 
handled  can  be  either  of  the  lumped  or  consistent  types. 
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The  algorithm  programmed  in  ADINA  is  the  determinant 
search  method  as  defined  in ^  ^  The  algorithm  combines  two 
basic  steps  namely 

i)  Triangular  factorization  and 
ii)  Vector  inverse  iteration. 

The  eigenvalues  are  obtained  in  a  sequence  starting  from  the 
least  dominant  eigenvalue  where  it  should  be  noted  that  the 
lowest  may  be  zero.  Because  of  this,  the  algorithm  can  be 
employed  to  calculate  the  rigid  body  modes. 

Once  one  eigenvalue/vector  is  obtained,  a  fairly  effici¬ 
ent  accelerated  secant  iteration  procedure  is  used  to  obtain  a 
shift  near  to  the  next  unknown  eigenvalue.  The  eigenvalue 
separation  theorem,  Sturm's  sequence  property,  is  employed 
in  this  iteration.  The  determinant  evalution  involves  a  tri¬ 
angular  factorization  of  the  eigenvalue  problem  defined  by 
Eq.  (II. 4).  Once  a  shift  to  the  neighborhood  of  unknown  eigen¬ 
value  is  obtained,  inverse  iteration  is  used  to  calculate  the 
eigenvector.  The  associated  eigenvalue  is  then  extracted  ac¬ 
curately  by  employing  the  Rayleigh  quotient  correction  to  the 
shift. 

The  output  options  available  to  the  user  are  fairly  exten¬ 
sive.  The  more  important  of  these  are  defined  by: 

1.  Any  number  of  frequencies  can  be  calculated  up  to 
number  of  system  freedoms  and 

2.  Cut  off  frequency  option. 
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Ill  Computational  Capabilities  and  Shortcomings 


To  determine  the  overall  capabilities  and  shortcomings 
of  the  static  and  dynamic  algorithmic  apparatus  inherent  to  ADINA, 
several  major  behavioral  traits  must  first  be  ascertained, 
namely : 

1)  The  generic  algorithmic  sensitivities; 

2)  The  overall  convergence  characteristics  and  lastly' 

3)  The  generic  types  of  solution  blow  ups  encountered. 

The  evaluation  of  such  algorithmic  characteristics  will  involve 
the  determination  of  the  effects  of  several  factors.  These 
include  such  items  as: 

i)  The  load  increment  size; 

ii)  The  time  step  size? 

iii)  The  material  properties; 

iv)  The  kinematic  nonlinearity,  etc. 

In  the  case  of  the  eigenvalue/vector  extraction  algorithms, 
the  major  characteristics  to  be  considered  include: 

1.  Convergence  properties; 

2.  Eigenvalue/vector  deterioration; 

3.  Multiplicity  and  separability. 

Since  the  major  thrust  of  ADINA  (1977)  is  its  capability 
to  handle  nonlinear  materials  and  large  deformations  kinematics, 
the  main  emphasis  of  the  algorithmic  check  out  will  involve  both 
the  static  and  dynamic  solution  capabilities.  The  approach  taken 
to  evaluate  the  properties  of  such  solution  branches  will  be 
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essentially  three  fold  in  particular: 
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1)  Establish  and  perform  bench  mark  tests;  2)  Perform 
parametric  studies  to  evaluate  algorithmic  sensitivities  and; 

3)  Establish  comprehensive  convergence  tests. 

As  will  be  seen  later,  the  main  purpose  of  the  bench 
marking  will  be  to  establish  the  algorithmic  sensitivities, 
convergence  characteristics  and  typical  solution  blow  ups 
generic  to  ADINA.  In  this  context,  the  benchmarks  will  be 
chosen  so  as  to  establish  the  characteristic  anomalous  solu¬ 
tion  behavior.  Such  behavior  will  be  accentuated  by  performing 
parametric  studies  involving  such  factors  as  load  step  size, 
time  step  size,  material  properties,  kinematic  nonlinearities, 
etc. 

The  main  purpose  of  the  comprehensive  convergence  testing 
will  be  to  monitor  solution  degradation  as  the  various  above 
mentioned  factors  are  parametrically  varied.  The  use  of  such 
convergence  tests  in  conjunction  with  the  benchmark  problems 
will  enable  the  characterization  of  the  overall  solution  patho¬ 
logies  inherent  to  ADINA. 

In  addition  to  the  foregoing  static,  dynamic  and  eigen¬ 
value  check  outs,  several  benchmarks  will  be  presented  which 
evaluate  the  pathological  characteristics  inherent  to  the  Gaus¬ 
sian  quadrature  of  the  isoparametric  element  stiffnesses.  This 
will  include  such  considerations  as  determining  how  the  degree 
of  accuracy  is  affected  by  the  integration  order. 

In  the  context  of  the  foregoing,  the  subsections  to  fol¬ 
low  will  include  discussions  on  the  establishment  of  comprehen¬ 
sive  convergence  tests,  as  well  as  on  the  results  of  parametric 
studies  involving  a  variety  of  benchmark  problems  including  static, 
dynamic,  eigenvalue  and  element  stiffness  integration  evaluations. 

III.l  Convergence  Tests 

To  establish  the  pathologies  of  the  algorithmic  sensi¬ 
tivities  of  ADINA,  a  comprehensive  convergence  test  which  can 
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monitor  both  the  quality  and  direction  of  a  given  iteration 
process  must  be  developed.  In  this  context,  the  standard  norm 
tests  on  the  nodal  increment  vector  and  out  of  balance  loads 
cannot  serve  in  such  a  capacity.  This  follows  from  the  fact 
that  the  basic  structure  of  such  tests  is  geared  primarily  to 
determining  outright  success  or  failure  of  solution  rather 
than  the  quality  of  convergence.  In  particular,  what  is  neces¬ 
sary,  is  the  ability  to  track  both  the  local  as  well  as  global 
behavior  of  the  iterated  solution. 

Before  developing  such  a  convergence  monitor,  it  is 
worthwhile  to  note  that  for  statically  loaded  problems,  the 
manner  of  loading  holds  the  key  to  the  direction  of  proper  con¬ 
vergence.  For  instance,  assuming  that  the  loading  is  mono- 
tonically  applied,  then  one  would  expect  essentially  monotone 
behavior  in  the  various  field  variables.  Similar  statements 
can  be  made  of  the  early  stages  of  the  transient  response  pro¬ 
blem.  Such  behavioral  traits  should  obviously  be  expected  of 
the  iteration  process.  Namely,  monotone  loading  should  always 
lead  to  monotone  decreasing  iterates  and  monotone  field  vari¬ 
ables.  As  it  is  computationally  inefficient  to  track  the 
monotonicity  of  the  various  field  variables  such  as  displace¬ 
ment,  stress  or  strain,  an  alternative  measure  should  be  em¬ 
ployed  . 

In  the  context  of  the  foregoing,  the  local  and  global 
strain  energy  functions  are  admirably  suited  for  the  current 
purpose.  This  follows  from  the  fact  that  in  addition  to  the 
all  important  global  and  load  monotonicity  properties,  the 
strain  energy  also  has  positive  or  negative  definite  attri¬ 
butes.  Namely,  as  the  calculations  proceed  for  a  monotone 


loading  situation,  one  would  expect  a  sequence  of  monotone, 
positive  definite  energy  iterates.  As  will  be  seen  in  the 
next  section,  behavior  to  the  contrary  will  almost  always  lead 
to  solution  divergence  or  at  best  poor  convergence  characte¬ 
ristics. 

Based  on  applications  of  the  virtual  work  principle  and 
the  incremental  approach,  the  strain  energy  stored  is  given 
by  the  expression  (Derivations  given  in  Appendix  Al) 
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where  Et  is  the  total  strain  energy  stored  during  an  itera¬ 
tion  cycle  involving  I  loadsteps,  k1  the  number  of  iterations 
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required  for  convergence  during  i  load  increment  and  AY^ 
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the  k  nodal  displacement  increment  of  the  i  load  step. 
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By  interpreting  the  nodal  deflection  vector  Y.  as  either  of 

^i 

an  element  or  global  form  (III.l)  can  be  used  to  describe 
the  local  element  or  global  strain  energies  stored.  Lastly, 
due  to  its  generalized  form  (III.l)  can  be  employed  for  all 
the  material  models  inherent  to  ADINA.  For  the  sake  of  con¬ 
venience,  the  derivation  of  the  expression  for  is  given  in 
Appendix  2. 

In  order  to  perform  energy  checks  of  the  solution, 
various  modification  had  to  be  implemented  in  the  ADINA 
program,  so  as  not  to  adversely  effect  core  requirements  and 
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dynamic  allocation,  a  new  direct  access  file  was  defined. 

To  increase  the  generality  of  the  energy  option,  the  addi¬ 
tion  was  programed  so  as  to  handle  the  following  possibilities: 

1)  Multiple  element  groups; 

2)  Multiple  solution  blocks; 

3)  All  material  property  groups  and; 

4)  All  algorithmic  options. 

III. 2  Benchmarking/Parametric  Studies 

While  the  benchmarking  discussed  in  this  section  serves 
the  purpose  of  a  check  on  the  accuracy  of  the  programming  and 
adequacy  of  the  theory  inherent  to  ADINA  (1977),  its  main 
function  is  to  determine  the  computational  capabilities  and 
shortcomings  of  the  algorithmic  apparatus.  For  convenience, 
those  benchmarks  which  were  run  to  strictly  check  out  the 
code  are  discussed  in  Appendix  2  and  the  input  data  are 
included  in  Appendix  A3.  In  this  context,  as  noted 

earlier,  this  section  will  overview  the  results  of  extensive 
parametric  studies  involving  the  effects  of  such  factors  as: 

i)  Load  step  size; 

ii)  Time  step  size; 

iii)  Material  properties; 

iv)  Kinematic  nonlinearities  etc. 

The  main  thrust  of  the  parametric  studies  will  be  to  determine 
the  pathology  of  anomalous  behavior  inherent  to  the  ADINA 
algorithmic  apparatus. 
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At  this  juncture,  it  should  be  noted  that  the  field 
theories  and  algorithmic  approaches  are  up-to-date  and  ac¬ 
curately  programmed.  Hence,  the  anomalous  behavior  depicted 
herein  rather  than  being  inherent  to  ADINA  are  generic  to 
nonlinear  structural  analysis  in  general. 

For  convenience,  the  overall  section  will  be  divided 
into  several  subsections  concentrating  on  such  items  as  the 
static,  dynamic  and  eigenvalue  branches  of  solution. 

III. 2a  Static  Solution 

This  subsection  will  consider  the  algorithmic  sensiti¬ 
vities  of  the  static  solution  branches.  Because  of  this  em¬ 
phasis,  the  geometric  configuration  of  the  benchmarks  con¬ 
sidered  will  purposely  be  kept  simple  yet  diverse  enough  to 
account  for  geometric  configuration  effects. 

Interestingly,  while  such  factors  as  geometry,  mate¬ 
rial  properties,  boundary  conditions,  etc.  all  have  some 
effect  on  the  choice  of  load  increment  size,  once  an  exces¬ 
sive  value  has  been  chosen,  typically  similar  types  of  solu¬ 
tion  degradation  are  encountered.  In  this  context,  the  first 
benchmark  chosen  will  be  used  to  establish  much  of  the  generic 
solution  anomalies  inherent  to  situations  wherein  excessive 
load  step  increments  are  considered.  Figure  III.l  illustrates 
the  geometry,  loading  and  material  properties  of  a  rubber 
sheet  which  serves  as  benchmark  1.  This  problem  combines  kine¬ 
matic,  kinetic  and  material  nonlinearity  since  large  defor¬ 
mations  and  Mooney  type  materials^  are  admitted  in  the  model. 
To  establish  the  generic  solution  anomalies  of  the  static 
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algorithm,  three  basic  scenarios  typical  of  ADINA  will  be 
excited  namely: 

1.  Outright  solution  failure  in  first  load  step; 

2.  Solution  failure  in  successive  load  steps; 

3.  Poorly  converged  solutions. 

The  actual  nature  of  solution  degradation  depends  on 
the  relative  magnitude  of  incrementation.  Typically  three 
basic  types  of  solution  pathology  occur.  These  can  be 
categorized  by: 

1.  Immediate  and  strong  nonmonotonicity ; 

2.  Moderate  but  progressively  increasing  non¬ 
monotonicity  and  nonpositive  definiteness  and; 

3.  Mild  monotonicity  with  either  very  gradual  in¬ 
creases  or  decreases  in  solution  oscillations. 

Note  such  behavior  can  be  excited  either  in  the  first  or 
successive  load  steps.  For  example.  Fig.  III. 2  illustrates 
immediate  nonmonotonic ity  in  the  first  load  step  caused  by 
excessive  incrementation  size.  While  solution  failure  for 
the  given  load  step  was  initiated  by  a  Fortran  overflow  stop, 
typically  out-of-balance  loads  are  encountered  for  such  pro¬ 
blems  . 

For  the  given  problem,  for  load  increments  somewhat 
lower  than  those  of  the  proceeding  example,  typically  mode¬ 
rate  but  progressively  increasing  nonmonotonicity  can  be 
excited  either  in  the  first  or  successive  load  steps.  Figu¬ 
res  III.  3-5  illustrates  such  solution  degradation.  As  can  be 
seen,  for  the  given  load  increment,  excellent  convergence  is 
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obtained  in  the  first  load  step.  In  the  second,  a  mild 
form  of  nonmonotonicity  and  nonpositive  definiteness  is  en¬ 
countered.  Since  this  increment  required  a  large  number 
iterations,  to  enable  the  solution  to  proceed  further,  the 
iteration  count  limit  was  raised.  In  this  context,  strong 
and  progressively  increasing  nonmonotonic ity  and  nonpositive 
definiteness  is  encountered  in  the  third  load  step  as  illu¬ 
strated  in  Fig.  III. 5.  Here  solution  failure  is  due  to  out-of¬ 
balance  loads. 

A  better  perspective  of  the  foregoing  behavior  can  be 
obtained  by  a  juxtaposition  of  results.  This  is  given  in  Fig. Ill 
6 .  The  scenario  of  solution  degradation  depicted  is  typical 
of  excessive  load  incrementation.  Note  as  can  be  seen  from 
these  results,  the  onset  of  such  behavior  is  signalled  by 
the  initiation  of  nonmonotonicity. 

In  the  preceding  discussion,  emphasis  was  given  to 
global  considerations.  Additional  insights  can  be  obtained 
by  studying  the  behavior  of  the  local  element  energies.  Such 
an  approach  will  enable  the  tracking  of  the  progressive  solu¬ 
tion  degradation  from  a  purely  local  point  of  view.  In  this 
context,  Figures  III. 7-13  illustrate  selected  local  conver¬ 
gence  characteristics  for  the  global  results  depicted  in 
Figures  III. 3-6.  As  can  be  seen,  the  solution  degradation  is 
initially  localized  but  gradually  spreads  to  the  entire  struc¬ 
ture  as  the  iteration  process  continues.  Because  of  this,  it 
appears  quite  plausible  that  the  choice  of  increment  size 
should  reflect  local  solution  characteristics  so  as  to  con¬ 
trol  the  growth  of  localized  degradation. 


Fig.  1 1 1  -  6  Global  Energy  Increment  of  Rubber 

Sheet  (1st,  2nd,  and  3rd  Load  Steps) 
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Figure  I II. 14  illustrates  the  geometry,  loading  and 
material  properties  of  a  fixed  arch  which  serves  as  bench¬ 
mark  2.  This  problem  involves  geometric  (kinematic)  non¬ 
linearity  which  is  strongly  influenced  by  the  geometry  of 
the  arch.  As  noted  earlier,  improper  incrementation  and 
algorithmic  choice  can  lead  to  solution  failure.  This  is 
particularly  true  for  arch  problems  which  are  subject  to 
snap  through  if  excessively  loaded.  While  anomalous  solu¬ 
tion  behavior  will  always  occur  in  the  vicinity  of  buckling 
and  bifurcation  points,  improper  incrementation  will  make  the 
accurate  estimation  of  such  points  impossible. 

As  with  the  previous  benchmark,  choosing  excessively 
large  load  steps  can  lead  to  various  types  of  solution  de¬ 
gradation.  For  the  current  example  problem,  such  behavior 
can  be  exacted  far  in  advance  of  bifurcation  or  buckling 
zones.  In  particular.  Figures  III. 15  and  III. 16  illustrate  typi¬ 
cal  types  of  solution  breakdowns  in  various  of  the  load  steps. 
Interestingly,  while  the  nature  of  the  nonlinearity  of  the 
arch  and  rubber  sheet  problems  are  inherently  different,  the 
pathology  of  solution  failure  appears  qualitatively  the  same. 

This  is  essentially  an  outgrowth  of  the  fact  that  the  MNR 
algorithm  is  slope  driven. 

In  addition  to  monitoring  global  energy,  local  ele¬ 
ment  characteristics  were  also  tracked.  Here  again,  solu¬ 
tion  degradation  was  found  to  be  initially  localized  but 
spreads  to  neighboring  elements  as  the  iteration  process 
continues.  The  rate  of  such  spreading  is  associated  with 
the  topology  of  the  solution  space.  In  this  context,  as  was 
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noted  earlier,  to  improve  local  solution  characteristics, 
the  load  incrementation  process  should  involve  the  use  of 
localized  information  concerning  element  stiffness  behavior 
namely  positive,  negative  definiteness, monotonicity ,  etc. 

Note,  waiting  until  local  solution  anomalies  spread  to  the 
entire  structure  may  lead  to  restart  difficulties.  Namely 
since  localized  solution  breakdowns  can  occur  earlier  in  the 
iteration  process,  all  successive  iterates  tend  to  drift 
from  the  true  solution.  Hence  upon  restart,  difficulties 
arise  as  to  what  portion  of  the  solution  is  valid.  Obviously, 
the  constant  monitoring  of  local  behavior  can  mitigate  such 
restart  problems. 

As  can  be  seen  from  the  foregoing  benchmarks,  for  static 
problems  excessive  incrementation  leads  most  frequently  to 
error  stops  involving  out-of -balance  loads  or  out  of  iteration 
counts.  Such  solution  difficulties  can  under  the  proper  circum¬ 
stances  be  circumvented  by  the  use  of  improved  iteration  pro¬ 
cedures.  For  example.  Figures  III. 17  and  18  illustrate  re¬ 
sults  obtained  for  benchmarks  1  and  2  wherein  an  NR  algorithm 
employing  constant  updating  was  used.  Note  while  the  modified 
NR  approach  yields  divergent  solutions  for  such  problems, 
the  constantly  updated  version  succeeded  in  either  proceeding 
further  or  outrightly  obtaining  the  final  solution.  An  al¬ 
ternative  to  the  straight  update  approach  is  yielded  by  al¬ 
lowing  the  solution  to  proceed  via  the  modified  INR  algorithm 
until  solution  degradation  is  initiated  either  locally  or 
globally.  Once  this  is  encountered,  the  algorithm  can  then 
be  switched  to  the  constant  update  approach  to  finish  out  the 
solution. 
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Further,  static  benchmarking  involving  different 
geometries  (beams,  plates,  spherical  shells/caps  and  mate¬ 
rial  properties  (plasticity)  also  initiated  similar  static 
solution  pathologies.  The  common  denominator  to  all  such 
solution  breakdowns  was  the  initiation  of  localized  non¬ 
monotonicity.  Once  encountered,  typically  such  behavior 
spreads  to  the  entire  structure  either  in  the  given  load 
step  or  in  successive  load  steps.  Note  such  solution  ano¬ 
malies  are  not  generic  to  ADINA  but  rather  are  an  outgrowth 
of  the  use  of  MNR  algorithms  which  are  based  on  tangent  stiff¬ 
ness  properties.  Hence,  all  general  purpose  nonlinear  codes 
employing  such  algorithmic  approaches  should  encounter  such 
anomalous  behavior  if  improper  load  incrementation  is  attempted. 
Section  IV  discusses  several  possible  improvements  which  can 
be  implemented  into  the  ADINA  program  to  diminish  such  dif¬ 
ficulties. 


III. 2b  Dynamic  Solution 

Before  considering  the  algorithmic  sensitivities  of 
the  dynamic  solution  branch,  several  important  points  must 
first  be  clarified.  The  most  important  of  these  is  the  need 
to  define  the  manner  in  which  the  inertial,  stiffness  and 
damping  effects  interact  to  control  the  overall  dynamic  res¬ 
ponse  of  a  structure.  As  a  starting  point,  we  shall  first 
consider  the  purely  linear  case.  Next,  the  effects  of  struc¬ 
tural  nonlinearity  (kinematic,  kinetic  &  material)  will  be 
considered.  Such  insights  will  enable  us  to  establish  the 
proper  framework  from  which  to  evaluate  the  dynamic  solution 
algorithm. 
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As  noted  earlier,  the  main  thrust  oi  this  work  will 
be  to  evaluate  the  algorithmic  sensitivities  of  ADINA  (1977) . 
As  with  the  static  branch,  this  will  be  accomplished  by 
exciting  the  various  forms  of  pathological  behavior  inherent 
to  the  algorithms. 

Since  the  dynamic  branch  of  ADINA  employs  direct  nu¬ 
merical  integration  to  yield  the  transient  solution,  the 
problem  of  central  importance  involves  the  matter  of  time 
step  size.  This  question  itself  consists  of  two  major  points 
namely: 

1.  System  dynamic  characteristics  and; 

2.  Characterization  of  exciting  fields. 

The  problem  of  dynamic  system  characterization  can  generally 
be  answered  by  considering  the  properties  of  its  natural 
frequencies.  This  includes  such  questions  as  spectral  spa¬ 
cing,  size  of  lowest  frequencies,  modal  participation  factors, 
etc.  In  the  case  of  linear  systems,  once  such  characteristics 
are  defined,  they  undergo  no  changes  during  the  calculation 
phased6  Vhis  is  obviously  a  direct  outgrowth  of  the  linearity 
of  the  governing  differential  operators  characterizing  such 
problems . 

The  characterization  of  the  exciting  fields  also  in 
essence  requires  the  definition  of  its  spectral  properties. 

In  this  context,  several  possible  situations  may  arise  namely: 

1.  A  profusion  of  frequencies  may  be  present; 

2.  Small  numbers  of  frequencies  are  present  or; 

3.  Selected  frequencies  dominate  behavior. 

Regardless  of  which  situation  exists,  once  the  spectral 


characteristics  of  the  structural  system  and  exciting  fields 
are  established,  the  proper  time  step  increment  can  be  sele¬ 
cted  so  that  the  participation  of  the  dominant  system  fre¬ 
quencies  can  be  accounted. 

For  situations  in  which  a  profusion  of  exciting  fre¬ 
quencies  exist,  the  explicit  procedure  must  be  employed.  In 
ADINA  the  central  difference  algorithm  can  be  used  to  handle 
such  situations.  In  the  case  where  small  numbers  of  fre¬ 
quencies  or  selected  frequencies  dominant  the  behavior,  im¬ 
plicit  operators  are  employed  to  handle  the  response.  For 
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such  situations,  the  Wilson  and  Newmark  integration  sche¬ 
mes  are  available  in  ADINA. 

Regardless  of  the  approach  employed,  once  the  full 
spectrum  of  excitation  and  system  frequencies  are  accounted 
for,  due  to  the  stationarity  of  the  spectral  properties  of 
linear  systems,  no  changes  need  be  implemented  in  the  time 
step  to  capture  the  response.  Such  is  not  the  case  for  non¬ 
linear  systems.  In  particular  for  nonlinear  systems,  major 
changes  in  the  spectral  properties  can  occur  from  moment  to 
moment  during  the  excitation.  Generally,  three  basic  types 
of  changes  can  occur  in  the  structural  behavior.  These  can 
be  categorized  by: 

i)  global/localized  structural  stiffening; 

ii)  global/localized  structural  softening  or; 

iii)  combined  structural  stif fening/softening . 

For  problems  which  are  characterized  by  material  nonlinearity, 
generally  softening  behavior  is  usually  encountered. 


A  good  example  of  such  behavior  can  be  seen  in  stru¬ 
cture  undergoing  small  plastic  type  deformations.  During 
the  development  of  plastic  flow,  such  structure  behave  softer 
than  their  purely  elastic  counterparts.  Interestingly,  while 
the  initial  time  step  chosen  for  the  elastic  version  of  such 
structure  can  capture  the  inertial  effects  of  either  the 
linear  or  nonlinear  situations,  it  is  possible  that  the  level 
of  nonlinearity  excited  by  plasticity  can  be  severe  enough 
to  cause  solution  difficulties.  Note  the  problems  which  are 
typically  encountered  for  such  situations  are  a  direct  out¬ 
growth  of  the  inability  to  handle  the  nonlinearity  via  MNR 
type  algorithms.  That  is,  if  extensive  softening  is  encoun¬ 
tered  during  a  given  time  step,  the  various  solution  patho¬ 
logies  depicted  for  the  static  case  can  be  excited.  This 
follows  from  the  fact  that  the  level  of  nonlinear  excited 
during  a  given  time  step  can  be  likened  to  that  encountered 
during  an  excessively  large  static  load  step.  In  either  case, 
the  solution  degradation  follows  the  same  pattern  as  the 
static  case  namely,  out  of  balance  loads  or  iteration  count 
stops  can  be  initiated.  Since  the  time  step  size  controls 
the  degree  of  nonlinearity  excited,  the  only  approach  avail¬ 
able  to  remedy  this  situation  is  to  decrease  the  time  step 
size.  In  ADINA  such  modifications  can  be  established  for 
selected  times  via  the  block  option.  If  the  adjustment  of 
time  step  size  cannot  stem  the  occurrence  of  out  of  balance 
loads  or  out  of  iteration  stops,  the  reformation  but  no  itera¬ 
tion  option  can  be  employed.  Solution  failure  under  this 
option  generally  involves  the  initiation  of  negative  pivots. 


This  follows  from  the  fact  that  without  the  out  of  balance 
load  check  inherent  to  the  iteration  option,  there  is  no 
bound  on  the  growth  of  load  imbalance. 

For  plastic  problems,  perhaps  the  most  severe  form  of 
material  nonlinearity  occurs  in  the  neighborhood  of  load 
unload  zones.  In  such  situations,  as  noted  in  Fig.  II. 8, 
transitions  between  softening  to  hardening  to  softening  be¬ 
havior  can  be  encountered.  The  major  solution  difficulty 
which  arises  from  such  behavior  is  the  fact  that  the  modified 
NR  algorithm  can  miss  the  need  to  unload  during  the  iteration 
phase  of  a  given  load  step.  This  usually  leads  to  solution 
drift  and  typically  out  of  balance  loads.  To  bypass  such 
difficulties,  the  only  avenue  open  to  the  ADINA  user  is  to 
stop  equilibrium  iteration  in  the  neighborhood  of  unloading. 
This  obviously  can  be  initiated  through  the  use  of  the  re¬ 
start  and  block  options.  While  such  an  approach  can  handle 
the  load/unload  problem,  the  apriori  knowledge  of  the  exis¬ 
tence  of  such  behavior  is  usually  not  immediately  available. 
Hence  such  information  must  be  ascertained  by  extensive  and 
oftentimes  costly  parametric  studies  involving  the  restart 
and  block  options. 

In  the  case  of  structure  with  severe  hardening  char¬ 
acteristics,  the  modified  NR  algorithm  inherent  to  ADINA  is 
destined  to  either  converge  poorly  or  initiate  out  of  balance 
loads  as  depicted  in  Fig.  III.  19.  The  only  way  to  avoid  such 
difficulties  is  to  employ  the  reformation  option  without  ite¬ 
ration.  This  obviously  can  lead  to  solution  drift  as  the 
solution  proceeds. 


Norm  of  Nodal  Displacements  1  1  -l  11 

Typical  Solution  Drift  for  Softening/Hardening 
Probl ems 
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In  the  context  of  the  foregoing  discussion,  the  bench¬ 
marking  of  the  dynamic  solution  branch  of  ADINA  discussed  in 
this  section  will  attempt  to  establish  the  various  forms  of 
pathological  behavior  intrinsic  to  the  direct  integration 
algorithm.  Note  the  main  emphasis  of  such  numerical  experi¬ 
ments  will  be  to  establish  the  important  effects  of  non¬ 
linearity  on  the  generation  of  anomalous  solution  behavior. 
This  emphasis  follows  from  the  fact  that  the  main  thrust  of 
ADINA  is  its  ability  to  handle  both  material  and  geometric 
nonlinearity.  As  with  the  static  branch,  because  of  this 
emphasis,  the  geometric  configuration  of  the  benchmarks  con¬ 
sidered  will  purposely  be  kept  simple  yet  diverse  enough  to 
account  for  geometric  configuration  effects. 

As  will  be  seen  from  the  following  discussion,  while 
such  factors  as  geometry,  material  properties  and  boundary 
conditions  all  have  varying  effects  on  the  choice  of  time 
step  size,  once  an  excessive  value  has  been  chosen,  typically 
similar  types  of  anomalous  solution  behavior  are  encountered. 
For  instance,  once  such  a  situation  is  excited,  the  iteration 
loops  tend  to  initiate  either  out  of  balance  load  or  out  of 
iteration  stops.  For  the  noniterative  loop,  typically  the 
negative  pivot  check  is  encountered. 

Figure  III. 20  illustrates  the  geometry  and  material 
properties  of  a  cantilevered  beam  which  serves  as  the  first 
benchmark  of  the  dynamic  solution  branch.  Since  large  de¬ 
formations  and  either  elastic  or  plastic  (strain  hardening) 
material  properties  are  employed,  this  problem  combines  the 
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effects  of  kinematic,  kinetic  and  material  nonlinearity. 

Since  we  are  primarily  concerned  with  the  characterization 
of  algorithmic  behavior,  to  establish  the  validity  of  model 
employed,  runs  were  performed  involving  20,  40,  80  and  160 
elements  for  selected  problems.  This  enabled  the  deter¬ 
mination  of  whether  modelling  difficulties  initiated  solution 
degradation  as  opposed  to  algorithmically  generated  ano¬ 
malous  behavior . 

To  establish  the  anomalous  behavior  of  the  various 
options  of  the  transient  algorithms,  the  problems  run  in¬ 
volved  increasing  orders  of  severity.  This  enabled  the 
evaluation  of  the  sources  of  degradation.  In  this  context. 
Figs.  III.  21-26  illustrate  the  various  characteristics  of 
the  displacements,  velocities  and  accelerations  associated 
with  the  large  deformation  responses  of  an  elastic  canti¬ 
levered  beam  subject  to  a  concentrated  step  load.  These 
were  run  with  the  Newmark  option.  Similar  results  were  ob¬ 
tained  via  the  Wilson  operator.  As  can  be  seen  from  the 
velocity  and  acceleration  fields  depicted,  because  of  the 
lack  of  damping,  numerous  higher  order  spectral  modes  par¬ 
ticipate  in  the  overall  response.  This  necessitates  the 

selection  of  an  extremely  small  time  increment  namely 

-4 

At'''.  5  x  10  sec.  While  such  behavior  eventually  initiates 
long  term  solution  drift,  for  the  first  few  cycles  of  time, 
good  accuracy  is  maintained. 

To  increase  the  degree  of  severity,  plasticity  is  ad¬ 
mitted  in  the  model.  As  noted  earlier,  this  initiates  sof¬ 
tening  type  behavior.  Namely,  the  apparent  period  of  the 
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response  tends  to  increase  because  of  the  softening  structural 
stiffness  initiated  by  plastic  transitions.  Since  such  a 
trend  tends  to  improve  the  inertia  following  characteristics 
of  the  numerical  integration,  this  enables  a  good  check  on 
the  stiffness  following  capabilities  of  the  algorithm.  The 
first  attempt  at  this  problem  employed  the  use  of  the  equi¬ 
librium  option  and  a  distributed  type  loading.  The  results 
are  depicted  in  Figs.  III. 27-32.  As  noted  earlier,  the  ite¬ 
ration  option  of  the  NR  algorithm  inherent  to  ADINA  does  not 
reform  the  stiffness  during  the  iteration  cycle.  Hence  upon 
unloading  which  typically  occurs  at  cycle  limits,  such  an 
algorithm  misses  the  transition  from  plastic  back  to  elastic 
during  the  iteration  phase.  This  can  lead  to  solution  de¬ 
gradation.  An  example  of  such  anomalous  behavior  is  il¬ 
lustrated  in  Figs.  III. 27-32.  The  actual  cause  of  failure 
is  the  initiation  of  the  out  of  balance  load  stop,  in  ADINA, 
such  difficulties  can  be  circumvented  by  employing  the  no 
iteration  option. 

Figures  III. 33-38  illustrate  the  elastic-plastic  (strain- 
harding)  distributed  load  response  as  obtained  by  the  no 
iteration  option.  As  can  be  seen,  the  imbalance  loads  encoun¬ 
tered  earlier  during  the  first  cycle  limit  were  not  monitored 
through  the  use  of  this  option.  Note  while  the  solution  ap¬ 
pears  successful  when  viewed  strictly  from  the  displacement 
point  of  view.  Figs.  III. 33  and  35,  such  a  simplistic  approach 
can  be  dangerous.  This  is  clearly  seen  by  considering  the 
behavior  of  the  velocity  and  acceleration  fields  depicted  in 
Figs.  III. 35-38.  As  can  be  seen  from  these  figures,  once 
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the  solution  passes  the  first  cycle  limit,  progressive  de¬ 
gradation  is  initiated.  This  follows  from  the  fact  that  with¬ 
out  iterative  corrections,  the  imbalance  loads  cause  a  drift 
on  the  solution  hyper  surface.  Because  of  this,  the  only  way 
the  solution  can  be  abnormally  terminated  in  ADINA  1977  is 
via  the  negative  pivot  stop.  Otherwise  the  calculations  will 
be  terminated  via  the  usual  time  limit  check. 

Increasing  the  problem  severity  further.  Figures  III. 

39-44  illustrate  the  elastic-plastic  cantilevered  beam  res¬ 
ponse  to  a  concentrated  load.  Due  to  the  severity  of  the 
loading,  outright  failure  occurred  when  the  iterative  solu¬ 
tion  option  was  employed.  To  circumvent  such  difficulties 
the  reformation  but  no  iteration  option  was  employed  to  gene¬ 
rate  the  results  depicted  in  these  figures.  Note  while  ade¬ 
quate  solution  behavior  is  obtained  during  the  early  stages 
of  the  response,  as  time  increases,  the  latter  phases  are 
marked  by  progressive  degradation.  This  can  be  seen  from  the 
growing  oscillations  which  occur  in  the  velocity  and  accele¬ 
ration  fields  shown  in  Figs.  III. 41-44.  These  culminate  in 
the  initiation  of  a  negative  pivot  stop.  As  noted  earlier, 
this  is  a  direct  outgrowth  of  the  fact  that  without  iterative 
corrections,  the  load  imbalance  tends  to  grow  as  the  time 
incrementation  proceeds. 

Figure  III. 45  illustrates  the  geometry  and  material  proper¬ 
ties  of  the  spherical  cap  which  serves  as  the  second  benchmark  of 
the  dynamic  solution  branch.  As  with  the  beam,  large  displace¬ 
ment  elastic-plastic  behavior  is  admitted  in  the  model.  Whereas 
the  previous  benchmark  established  solution  difficulties 
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associated  with  long  term  behavior  (several  cycles),  the 
current  problem  will  emphasize  anomalous  behavior  during  the 
initial  stages  (first  cycle)  of  solution. 

In  this  context,  the  effects  of  severe  nonlinearity  on 
the  choice  of  time  step  size  will  be  considered.  To  achieve 
this  end,  the  spherical  cap  will  be  excited  by  step  pressure 
loadings  of  increasing  magnitude.  The  various  effects  of 
kinematic/material  nonlinearity  on  solution  success  is  shown 
in  Fig.  III. 46-49.  Note  as  illustrated,  three  basic  stages  are 
seen  to  exist.  These  can  be  categorized  by: 

1)  "Mild"  nonlinear  effect; 

2)  "Moderate"  nonlinear  effects  and; 

3)  "Strong"  nonlinear  effects. 

These  stages  of  behavior  will  be  considered  in  the  context  of 
time  step  sizes  which  are  successful  in  describing  the  linear 
response.  Before  discussing  such  behavior,  it  should  be  noted 
that  for  loads  below  the  critical  buckling  value,  the  cap 
exhibits  both  kinematic  and  material  softening. 

As  noted  earlier,  time  step  sizes  adequate  to  capture 
the  inertial  effects  of  linear  problems  will  remain  so  for 
softening  situations.  For  such  nonlinearities,  degradation 
typically  arises  as  an  outgrowth  of  the  inability  to  track 
the  changes  in  structural  stiffness.  For  mildly  nonlinear 
situations,  since  only  minor  changes  occur  in  the  structural 
stiffness,  typically  the  linear  time  step  is  adequate.  This 
can  be  seen  from  the  results  illustrated  in  Fig.  III. 46. 

Employing  the  same  time  step  size,  Fig.  III. 47  illustrates 
the  transient  cap  response  for  moderate  loadings.  For  such 
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situations,  the  higher  order  spectral  modes  and  the  material/ 
geometric  nonlinearity  combine  to  cause  an  extensive  amount 
of  iterating  to  capture  the  overall  solution.  Because  of 
this,  as  can  be  seen  from  Fig.  47,  solution  failure  is  ini¬ 
tiated  because  of  out  of  iteration  stops.  For  such  situations, 
generally  any  attempt  to  increase  the  number  of  iterations 
tends  to  fail  since  typically  an  extremely  large  number  of 
steps  is  required.  This  is  partially  an  outgrowth  of  the  local 
wavyness  of  the  solution  space  caused  by  the  presence  of  higher 
order  spectral  modes. 

For  the  given  time  size,  a  further  increase  in  loading 
tends  to  excite  significant  amounts  of  plastic  flow.  Since 
such  behavior  causes  extensive  amounts  of  structural  soften¬ 
ing  and  irreversible  energy  conversion,  much  of  the  higher 
order  spectral  modes  tend  to  be  washed  out  before  unloading 
occurs  at  the  first  cycle  limit.  Because  of  this,  interest¬ 
ingly  the  solution  tends  to  progress  further  before  the  out 
of  iteration  stop  is  excited.  This  behavior  is  clearly  il¬ 
lustrated  in  Fig.  48. 

The  three  stages  of  solution  noted  earlier  are  clearly 
seen  from  the  juxtaposed  results  depicted  in  Fig.  49.  Such 
behavior  was  obtained  for  both  the  Newmark  and  Wilson  inte¬ 
gration  options. 

Note,  further  dynamic  benchmarking  involving  different 
geometries  and  material  properties  tended  to  initiate  similar 
solution  pathologies  to  those  depicted  by  the  foregoing  pro¬ 
blems.  As  with  the  static  solution  branch,  such  anomalous 
behavior  does  not  appear  to  be  generic  to  ADINA.  Rather  it 
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appears  to  be  intrinsic  to  the  tangent  stiffness  algo¬ 
rithm  itself.  In  particular,  this  algorithm  appears  to  be 
sensitized  by  the  interactions  between  structural  nonlinearity 
and  inertial  effects.  For  moderately  nonlinear  problems,  such 
factors  combine  cause  a  "wavy"  solution  space  which  tends  to 
sensitize  the  tangent  stiffness  approach.  For  materials 
which  are  highly  dissipative,  such  behavior  is  somewhat  smo¬ 
othed  out  as  the  nonlinearity  is  further  increased. 

III. 2c  Eigenvalue/Vector  Extraction  Algorithms 

The  main  thrust  of  this  subsection  is  to  determine  the 
computational  characteristics  of  the  eigenvalue/ vector  ex¬ 
traction  algorithms  inherent  to  ADINA  (1977) .  In  particular 
this  will  include  considering  such  factors  as  the  convergence 
properties,  eigenvalue/vector  deterioration  and  multiplicity/ 
separability  characteristics. 

To  achieve  these  objectives,  the  benchmarking  employed 
must  be  capable  of  testing  the  extraction  algorithms  for 
situations  involving  arbitrary  eigenvalue  spacing.  In  this 
context,  the  algorithmic  check  out  considered  the  following 
types  of  situations  namely 

1.  Widely-spaced  frequencies; 

2.  Closely-spaced  frequencies; 

3.  Widely  space  groups  of  closely  spaced  frequencies 
and; 

4.  Multiple  frequency  branches 

The  foregoing  scenario  of  frequency  spacings  are  typical 
of  modern  day  support  structure.  The  most  complicated  of 
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these  usually  involve  main  structures  which  have  soft  or 
stiff  substructural  components.  Typically  such  components 
tend  to  lead  to  the  occurrence  of  widely  spaced  groups  of 
closely  arranged  frequency  eigenvalues.  For  instance,  sup¬ 
ported  symmetric  structures  are  good  examples  of  such  a 
situation.  This  follows  from  the  fact  that  symmetric  struc¬ 
ture  in  themselves  usually  have  multiple  frequency  branches. 
Supporting  such  structure  on  soft  foundation  supports  tends 
to  lead  to  an  overall  frequency  spectrum  involving  essent¬ 
ially  two  groups  of  eigenvalues,  namely; 

1.  Pseudo  rigid  body  modes  involving  the  main 
structure  acting  rather  rigidly  on  soft  supports 
and; 

2.  Main  structural  modes  which  have  been  separated 
by  the  presence  of  the  supports. 

Since  such  structure  can  be  made  to  exhibit  all  of  the 
eigenspacinqs  denoted  by  1)  -  4) ,  the  family  of  benchmarks 
employed  herein  involved  such  a  type  structural  arrangement. 

Before  describing  the  actual  benchmarks  employed,  it 
is  appropriate  to  note  that  in  order  to  quantify  potential 
eigenvalue/vector  deterioration,  a  fairly  extensive  number 
of  dependable  eigenvalues  must  be  available.  In  this  con¬ 
text,  the  approach  taken  herein  was  to  chose  a  structure 
which  could  be  solved  analytically  and  hence  have  any  number 
of  known  eigenvalues  and  associated  eigenvectors. 

Figures  III. 50-52,  show  the  overall  family  of  con¬ 
figurations  used  for  the  evaluation  of  the  eigenvalue/vector 
extraction  algorithm.  As  can  be  seen,  the  main  structure  of 
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v=.3  h  =  .  1  i  n  c  h 
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Fi g . 1 1 1 . 51  Geometry  of  Symmetrically  Supported 
Ring  on  Founda t i on 
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the  family  consists  of  a  ring  or  long  cylinder.  Fig.  III. 50. 

By  attaching  it  elastically  to  ground  via  a  series  of  radial 
or  tangential  springs,  a  variety  of  symmetric  or  asymmetric 
structure  can  be  generated,  Figs.  III. 51  and  52. 

Note  since  the  cylinder  is  itself  axisymmetric,  it 
thus  has  multiple  frequency  eigenvalues  of  multiplicity  two. 
The  spacings  of  the  multiple  pairs  of  eigenvalues  of  this 
structure  can  be  made  to  be  either  closely  or  widely  spaced 
depending  on  its  radius  to  thickness  ratio.  Such  spacings 
can  be  enhanced  by  attaching  the  cylinder  to  ground  via  a 
symmetric  radial  spring  support.  Fig.  III. 51. 

In  terms  of  this  structure,  the  problem  of  widely  spaced 
groups  of  closely /widely  spaced  frequency  pairs  can  be  gene¬ 
rated  by  introducing  asymmetrically  placed  tangential  spring 
supports  as  depicted  in  Fig.  III. 52.  Through  the  introduc¬ 
tion  of  the  asymmetric  support  springs,  the  multiple  frequency 
pairings  can  be  made  to  bifurcate  into  separate  frequency 
branches.  The  spacing  between  such  branches  can  be  made  to 
increase  or  decrease  depending  on  the  amount  of  asymmetric 
stiffness  introduced. 

The  choice  of  the  foregoing  benchmark  family  enabled  a 
logical  and  organized  approach  to  be  taken  for  the  problem 
of  evaluating  convergence  properties,  eigenvalue/vector  de¬ 
terioration  and  multiplicity/separability  characteristics. 

In  performing  the  aforementioned  benchmarking  activities, 
the  following  factors  were  monitored  during  the  calculation 
flow,  namely: 
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1.  The  required  number  of  iterations  per  eigenvalue; 

2.  The  convergence  characteristics  of  the  eigenvalues 
and  eigenvectors  and; 

3.  The  overall  solution  pathologies  as  the  iteration 
process  proceeded  from  low  to  high  frequencies. 

The  first  benchmark  employed  to  study  the  capabilities 
of  the  extraction  algorithm  consisted  of  the  cylinder  on 
radial  foundation  depicted  in  Fig.  III. 51.  In  addition  to 
monitoring  convergence  and  deterioration  properties,  this 
benchmark  enabled  the  evaluation  of  the  multiplicity  char¬ 
acteristics.  To  determine  the  eigenvalue/vector  deterio¬ 
ration  pathology,  the  cylinder  was  modelled  by  various  number 
of  elements  involving  either  80,  160,  300,  600  or  1000  nodes. 
Figure  III. 53  illustrates  a  typical  mesh  employed  in  the 
parametric  study.  Based  on  such  models,  the  eigenvalues  and 
associated  vectors  were  calculated.  In  all  the  cases  studied, 
the  convergence  characteristics  were  excellent.  Additionally, 
once  eigenvalue  or  vector  deterioration  was  noted  by  com¬ 
parisons  to  theoretical  values,  such  pathologies  could  almost 
always  be  directly  attributed  to  mesh  inadequacy  for  the  given 
mode . 

For  example,  natural  frequencies  were  obtained  for  a  160-element 

model  (Fig.  III. 53)  to  illustrate  the  deteriorated  multiple  root  eigen¬ 
value  family.  The  mode  shapes  associated  with  these  multiple 

eigenvalues  are  depicted  in  Figs.  III. 54-67.  As  can  be  seen 
from  the  results  illustrated,  the  modal  and  eigenvalue  dete¬ 
rioration  can  be  directly  attributed  to  the  lack  of  mesh  ade¬ 
quacy  for  higher  order  modes. 


The  only  difficulties  encountered  with  this  benchmark 
involved  problems  with  large  numbers  of  elements.  For  such 
cases,  as  noted  by  the  ADINA  manual  itself,  the  out  of  core 
solution  proved  somewhat  costly.  This  can  be  directly  at¬ 
tributed  to  the  necessity  of  using  low  speed  storage  (disks) 
reads  and  writes  to  pack  and  unpack  the  various  blocks  of 
the  global  stiffness  matrix. 

Figure  III. 52  illustrates  the  asymmetrically  supported 
cylinder.  By  adjusting  the  various  foundations  of  this 
model,  the  eigenvalue  spacing  was  varied  over  an  extensive 
range  of  separations.  This  enabled  a  further  evaluation  of 
the  multiplicity  and  separability  characteristics.  As  in 
the  previous  case,  the  cylinder  was  modelled  by  various  num¬ 
bers  of  elements  involving  either  80,  160,  300,  600  or  1000 
nodes.  In  the  parametric  studies  performed  with  these  models, 
the  convergence  characteristics  were  good  for  all  the  ranges 
of  separability  tested.  Whenever  eigenvalue  or  vector  dete¬ 
rioration  was  encountered,  it  could  generally  be  attributed 
to  the  lack  of  model  accuracy.  Namely,  the  lack  of  proper 
numbers  of  elements  required  to  adequately  define  a  given 
mode  shape . 

As  a  consequence  of  the  foregoing  parametric  benchmark 
studies,  it  follows  that  the  ADINA  eigenvalue/vector  extrac¬ 
tion  algorithm  has  good  convergence  characteristics  for  most 
typically  occurring  ranges  of  frequency  spacing.  The  only 
difficulties  encountered  were  associated  with  out  of  core 
solutions  wherein  slow  secondary  storage  was  involved.  Note 
by  increasing  the  size  of  the  dynamic  allocation  array,  larger 
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in-core  problems  can  be  handled  thus  enabling  raster  run 
times.  This  raises  an  interesting  point  concerning  the 
ADINA  overlay  structure  and  the  MVS  environment  typical 
of  IBM  systems.  While  such  studies  were  not  run,  it  would 
be  interesting  to  determine  the  advantages  of  totally  in 
core  Vrs  out  of  core  solutions  in  an  MVS  storage  environment. 
In  this  regard,  it  was  found  that  for  the  benchmarks  run 
for  the  study,  the  purely  in  core  approach  involving  large 
storage  dynamic  allocations  ran  faster  than  small  dynamic 
allocations  involving  the  out  of  core  option.  While  this 
might  seem  obvious  at  first,  it  must  be  remembered  that  in 
the  MVS  environment,  portions  of  the  program  are  paged  in  and 
out  of  the  central  core  during  the  course  of  running.  Such 
an  operation  can  be  likened  to  a  machine  generated  overlay 
structure.  In  this  context,  it  would  be  worthwhile  to  de¬ 
termine  in  future  studies  the  pros  and  cons  of  in  core  vrs 
out-of-core  program  architecture  in  an  MVS  environment. 


III. 2d  Numerical  Integration  of  Element  Stiffness 


In  the  evaluation  of  element  stiffness  for  isoparametric 
elements,  the  general  practice  is  to  use  Gauss  quadrature  nu¬ 
merical  integration.  Ordinarily,  the  degree  of  accuracy  in 
numerical  integration  is  affected  by  the  integration  order  cho¬ 
sen  consistent  with  the  approximating  shape  functions  used  for 

the  element.  For  example,  for  an  element  of  quadratic  displace- 

[71 

ment  field,  the  exact  numerical  integration  order  is  "3" 

However,  the  choice  of  integration  order  is  also  complicated  by 
other  factors,  i.e.,  variation  of  material  properties  over  the 
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element  (elastic-plastic  material) ,  incompressibility  constraints, 
presence  of  shear  stresses  in  thin-wall  structure,  etc.  In  this 
section,  some  numerical  difficulty  of  the  3/D  continuum  element 
is  revealed  when  it  is  used  for  bending  analysis. 

The  problem  considered  is  a  square  plate  simply  supported 
along  its  four  edges,  and  constrained  from  horizontal  motion  in 
the  directions  perpendicular  to  the  edges  of  the  supports.  The 
plate  is  subjected  to  a  uniformly  distributed  load.  Material 
properties  were  assumed  to  be  linearly  elastic  and  isotropic,  i.e. 

Young's  Modulus  E  =  2.1  x  106  Kg/cm2 

Poisson  ratio  v  =  0.3 

The  20  node  3/D  continuum  element  was  used  to  model  the  bending 
action  of  the  plate.  From  symmetry,  only  one  quarter  of  the 
plate  was  modeled  by  9  elements  with  96  nodes  as  shown  in  Fig. 

68.  The  plate  was  loaded  well  into  the  large  deflection  range 
so  that  membrane  action  of  the  plate  became  predominating.  Seve¬ 
ral  computer  runs  were  made  by  varying  the  thickness-ratio  (or 
aspect  ratio,  t/a)  of  the  plate  ranging  from  0.1  to  10  5.  For 
each  thickness  ratio,  two  numerical  integration  orders  were 
chosen  for  the  evaluation  of  element  stiffness;  namely,  3x3x2 
(exact)  and  2x2x2  (reduced)  orders,  two  integration  points  in 
the  thickness  direction  for  both  cases.  The  load  -  deflection 
responses  of  the  plate  with  two  different  aspect  ratios  are  plotted 
in  Figs.  69  and  70.  In  the  analysis,  the  results  obtained  from 
the  exact  integration  order  represent  convergent  solution  under 
equilibrium  iterations  are  totally  erroneous  as  seen  in  Figs. 

69  and  70.  On  the  other  hand,  the  reduced  integration  scheme 


gave  a  more  satisfactory  result  as  compared  to  the  analytical 

r  o  i 

solution  J •  The  non-dimensionalized  maximum  deflection  of 

Pa  4 

the  plate  vs.  various  aspect  ratios  for  the  load  factor  ^£-^-4.6,50 
and  208  are  plotted  in  Figj.  71  and  73,  respectively.  From  these 


plots,  one  can  clearly  see  that  the  use  of  3/D  continuum  elements 
with  exact  integration  order  fails  to  give  any  reasonable  solu¬ 
tion  for  thin  plates  with  a  n  aspect  ratio  smaller  than  10  3 . 

In  order  to  overcome  this  numerical  problem,  reduced  integration 
order  must  be  used  for  either  small  deformation  or  large  defor¬ 
mation  analysis.  It  is  also  noted  that  as  the  plate  thickness 

—  4 

becomes  very  thin  (t/a  less  than  10  ) ,  no  solution  can  be  ob¬ 

tained  even  the  reduced  integration  was  used. 

IV.  Potential  Algorithmic  Improvement 

Finite  element  simulations  of  statically  and  dynamically 
loaded  structures  composed  of  general  materials  undergoing 
large  displacements  usually  lead  to  nonlinear  field  equa¬ 
tions.  Since  the  types  of  nonlinearity  exhibited  by  such 
field  equations  are  both  diverse  and  complex,  the  question 
of  the  best  choice  of  an  appropriate  solution  algorithm  in¬ 
evitably  arises.  While  many  alternatives  are  available,  gene¬ 
rally  the  various  solution  procedures  may  have  special  ad¬ 
vantages  for  certain  classes  of  problems  but  may  exhibit  poor 
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convergence  for  other  situations.  For  example,  the  ideal 
general  purpose  nonlinear  finite  element  code  should  have 
several  algorithmic  options  augmented  with  a  degree  of  arti¬ 
ficial  intelligence.  Namely,  the  problem  solving  capability 
should  involve  a  heuristically  guided  trial  and  error  search 
in  the  space  of  possible  solution  via  an  automatically  struc¬ 
tured  algorithm.  Unfortunately,  currently  available  general 
purpose  finite  element  codes  of  which  ADINA  is  an  example,  do 
not  have  heuristic  (self-adaptive)  capabilities.  In  parti¬ 
cular,  because  of  its  wide  applicability,  like  ADINA  most 
such  codes  employ  some  primitive  (non-self-adaptive)  variant 
of  either  the  full  or  modified  incremental  Newton  Raphson 
(MNR)  algorithmic  procedure.  This  applies  to  both  the  static 
and  transient  solution  branches.  In  particular,  the  various 
transient  algorithms  such  as  Wilson,  Newmark,  Houbolt  and 
central  difference  are  usually  developed  around  the  primitive 
MNR  algorithm  and  hence  do  not  in  themselves  possess  heuristic 
(self-adaptive)  capabilities. 

While  the  ADINA  code  presents  the  user  with  far-reaching 
and  accurately  programed  capabilities,  without  apriori  physical 
insight,  expensive  parametric  studies  are  oftentimes  necessary  to 
assure  adequate  solution  convergence.  This  is  particularly  true 
for  unqualified  users  who  do  not  have  the  proper  training  in  non¬ 
linear  mechanics  theory.  As  has  been  seen  in  the  previous  sect¬ 
ion,  unless  the  proper  load  or  time  incrementation  is  employed, 
either  poor  convergence  or  out-of-balance  loads  are  generally  en¬ 
countered.  Such  anomalous  solution  behavior  is  not  an  out¬ 
growth  of  improper  programing  practice  or  illchosen  field 


representation,  but  rather  is  generic  to  all  nonlinear  F.E. 
codes  employing  non-self-adaptive  MNR  algorithms.  In  view 
of  such  shortcomings,  this  section  will  overview  a  series  of 
suggested  changes  to  the  ADINA  algorithmic  strategies  aimed 
at  increasing  the  efficiency  of  the  code  as  well  as  decreasing 
the  need  for  direct  user  intervention  which  is  generally 
quite  time-consuming  and  costly. 

From  an  overview  point  of  view,  the  main  thrust  of  the 
discussion  is  to  establish  a  three  level  strategy  in  parti¬ 
cular  : 


i) 

Level 

1: 

Preliminary  solution  development  via 
primitive  algorithms; 

ii) 

Level 

2: 

Solution  monitoring  via  validity/con- 

vergence  tests  and; 

iii) 

Level 

3: 

Self-adaptive  strategies  (Heuristic 
programming) . 

Like  the  current  ADINA  equation  solver,  the  first  level  should 
employ  primitive  (  MNR,  NR  etc.)  operators  to  generate  the 
solution  in  the  usual  manner.  The  second  level  should  involve 
the  constant  monitoring  of  the  different  stages  of  solution 
via  a  variety  of  validity /convergence  tests.  The  last  level 
should  consist  of  various  self-adaptive  strategies  which  are 
triggered  by  the  tests  initiated  by  the  level  two  surveillance 
process.  Note  these  strategies  should  be  arrayed  so  that  the 
correct  validity /convergence  tests  initiate  the  strategy  with 
the  appropriate  degree  of  sophistication. 

To  quantify  the  particulars  of  the  suggested  algorithmic 


improvements,  the  following  subsection  will  briefly  outline 
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the  main  intent  of  the  various  phases  of  solution  strategy 
noted  above. 

IV. 1  Level  1:  Primitive  Operators 

As  noted  earlier,  the  first  level  should  employ  pri¬ 
mitive  operators  to  generate  the  solution.  So  as  not  to 
cause  major  changes  in  the  ADINA  architecture,  such  operators 
should  involve  various  versions  of  the  NR  algorithm  namely: 

1.  Full  NR  with  constant  reformation  of  tangent 
stiffness  matrix  during  iteration; 

2.  nr  algorithm  with  intermittent  updating  during 
iteration; 

3.  Modified  NR  algorithm  with  tangent  stiffness  re¬ 
formation  at  beginning  of  load  step  but  not  during 
iteration  (currently  available) ; 

4.  NR  algorithm  with  iteration  but  no  reformation 
(currently  available) ; 

5.  NR  algorithm  with  reformation  but  no  iteration 
(currently  available)  and  lastly; 

6.  NR  algorithm  with  no  reformation  and  iteration 
(currently  available) . 

In  the  context  of  Fig.  II. 1,  Fig.  IV. 1  defines  the 
overall  logic  needed  to  implement  the  calculation  flow  as¬ 
sociated  with  the  NR  adaptions  noted  by  items  1-6  listed 
above . 

Such  an  addition  could  be  achieved  through  minor  modi¬ 
fication  to  the  EQUIT  subroutine  wherein  the  iterative  loop 
is  programmed.  Additionally,  new  control  parameters  would  be 
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necessary  in  such  subroutines  as  TDFE ,  QUADS  etc.  to  allow 
for  constant  and/or  intermittent  updating  of  the  tangent 
stiffness  matrix.  As  most  of  all  the  necessary  arrays  of 
data  are  available  in  EQUIT,  no  real  changes  in  code  archi¬ 
tecture  would  be  involved  apart  from  the  incorporation  of: 

1.  New  variables  to  control  the  choice  of  the  appro¬ 
priate  NR  algorithm; 

2.  New  logic  to  allow  for  the  proper  packing  of  the 
reformed  stiffness  in  secondary  storage. 

IV .  2  Level  2:  Validitv/Convergence  Tests 

The  validity/convergence  tests  should  be  the  core  of 
the  modification  to  ADINA.  For  the  present  purposes,  the 
discussion  will  be  organized  in  three  main  categories  namely: 

1.  Classical  convergence  tests; 

2.  "Quality"  of  convergence  tests; 

3.  Degree  of  "nonlinearity  tests". 

The  first  group  should  be  of  the  classical  normed  type 
pass  or  fail  variety  as  typified  by  the  convergence  test 
currently  programmed  in  ADINA,  namely: 

i)  the  out  of  balance  norm  test; 
ii)  nodal  displacement  norm  test. 

The  main  intent  of  such  tests  is  essentially  to  monitor  the 
success  or  failure  of  the  iterative  procedure.  While  such 
tests  are  efficient  and  well  adapted  to  this  prupose,  they 
cannot  really  forecast  potential  difficulties  until  outright 
failure  occurs.  Because  of  this,  if  the  restart  option  is 
employed,  depending  on  the  problem,  much  of  the  solution 
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generated  up  to  the  failure  point  may  have  to  be  scratched 
so  as  to  proceed  further.  In  this  context,  what  is  necessary 
are  so-called  validity  checks  which  enable  a  constant  moni¬ 
toring  of  the  solution  so  as  to  determine  whether  the  direc¬ 
tion  of  convergence  is  proper. 

In  the  context  of  the  foregoing,  the  second  group  of 
tests  should  be  concerned  with  the  quality  of  convergence 
namely  the  rate,  monotonicity,  positive,  negative  and  semi¬ 
definiteness,  etc.  Such  information  should  obviously  be  used 
to  trigger  various  modifications  in  the  primitive  iterative 
strategy  depicted  in  level  1.  As  can  be  seen  from  the  re¬ 
sults  outlined  in  Section  III,  various  aspects  of  the  system 
energy  could  be  employed  to  serve  as  validity  checks. 

In  addition  to  serving  as  a  quality  of  convergence 
test,  monitoring  various  aspects  of  the  energy  should  serve 
as  a  good  measure  of  the  degree  of  nonlinearity  excited  by 
successive  load/time  increments.  The  importance  of  such  tests 
follows  from  the  fact  that  although  finite  element  simula¬ 
tions  of  structures  composed  of  general  media  undergoing  large 
deflections  are  inherently  nonlinear,  the  degree  of  non¬ 
linearity  excited  varies  from  point  to  point  as  well  as  from 
load/time  increment  to  load/time  increment.  As  it  is  pos¬ 
sible  that  large  portions  of  a  structure  may  exhibit  basically 
linear  behavior,  the  third  phase  of  quality/ convergence  test¬ 
ing  will  enable  the  automatic  partitioning  of  the  structure  by 
allowing  for  preferential  updates  of  the  tangent  stiffness  de¬ 
pending  on  the  amount  of  local  nonlinearity  excited.  Such 
preferential  updating  obviously  represents  a  more  significant 
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modification  in  ADINA' s  architecture  than  the  first  two 
phases  of  quality /convergence  testing.  In  particular,  to 
accomodate  the  preferentially  updating  of  particular  stiff¬ 
ness  partitions,  such  a  change  would  entail  a  modification 
in  the  dynamic  allocation  scheme. 

IV. 3  Level  3:  Heuristic/Self-adaptive  Strategies 

Since  the  core  of  the  ADINA  solution  strategy  is  based 
on  the  tangent  modulus  approach,  any  Heuristic  programming 
modifications  should  obviously  be  centered  around  the  NR 
family  of  algorithms.  In  the  context  of  the  inherent  fea¬ 
tures  of  the  NR  algorithm,  the  adaptive  strategy  should  in¬ 
corporate  one  or  more  of  the  following  basic  options  namely: 

1)  Adaptive  tangent  stiffness  updates; 

2)  Adaptive  incremental  load  adjustments; 

3)  Adaptive  incremental  time  adjustments. 

As  noted  earlier,  the  various  self  adaptive  strategies 
should  be  programmed  so  as  to  be  triggered  by  various  of  the 
previously  discussed  validity/convergence  tests.  The  hie¬ 
rarchy  of  such  modifications  should  be  arrayed  so  that  the 
quality  of  convergence  test  can  be  used  to  initiate  the  al¬ 
gorithm  with  the  requisite  degree  of  sophistication  and  ef¬ 
ficiency  required  by  the  solution  difficulties  encountered. 

The  adaptive  stiffness  updates  can  themselves  be  or¬ 
ganized  in  three  main  categories  namely: 

a)  Global  updates, 

b)  Local  updates  and 

c)  Partial  updates. 


_ _ _ V _  _  .  - ...  ■  _ ..  ■fair  it  i 
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The  global  stiffness  updates  should  be  based  on  either  the 

degree  of  global  geometric  and  material  nonlinearity  excited 

or  because  of  poor  solution  convergence  as  determined  by  the 

validity /convergence  test.  The  initiation  of  local  updates 

should  depend  on  the  degree  of  local  element  nonlinearity 

excited.  The  incorporation  of  the  local  update  option  would 

obviously  enhance  the  overall  efficiency  of  the  calculations. 

While  the  global  and  local  preferential  updates  represent 

higher  order  solution  fixes,  pseudo  updates  such  as  the  BFGS 
r  9  1  * 

may  represent  a  less  costly  self-adaption  for  certain 
classes  of  problems  wherein  only  minor  changes  in  stiffness 
may  be  encountered.  Obviously  for  such  problems  as  plastic 
loading  and  unloading,  pseudo  updates  such  as  the  BFGS  are 
apt  to  yield  erroneous  results. 

The  adaptive  incremental  load  and  time  adjustments  can 
also  be  organized  in  several  main  categories  namely: 

a)  Increment  expansion; 

b)  Increment  contraction; 

c)  Corrective  incrementation. 

Increment  expansion  for  either  static  or  dynamic  problems 
can  be  initialized  by  various  of  the  quality  of  convergence 
checks.  Since  the  solution  quality  can  be  constantly  moni¬ 
tored,  any  degradation  initiated  by  increment  expansion  can 
be  followed  by  a  subsequent  contraction.  The  contraction  pro¬ 
cess  should  iteself  be  based  on  various  checks  such  as: 

1)  Poor  quality  of  convergence  (non-monotonicity; 
nonpositive  definiteness,  etc.); 

♦Currently  available  in  ADINA  1978 


2)  Level  of  nonlinearity  excited; 

3)  Nearness  to  bifurcation  points  as  monitored  by 
nonlinearity  and  energy  checks. 

Note,  a  high  level  of  local  nonlinearity  may  be  used  to 
initiate  local  modification  of  incrementation.  Obviously 
in  order  to  do  this,  special  bookkeeping  routines  must  be 
initiated  to  keep  track  of  the  status  of  the  loading. 

Note  for  load  increments  for  which  significant  solution 
divergence  is  encountered,  corrective  load  incrementation 
should  be  initiated.  In  particular,  negative  load  incre¬ 
mentation  should  be  employed  to  retrace  a  portion  of  load 
history  wherein  a  lower  order  algorithmic  strategy  yielded 
poorly  converged  results.  For  time  dependent  problems,  cor¬ 
rective  incrementation  should  be  initiated  to  retrace  a  por¬ 
tion  of  time  history. 


V. 


Summary 


As  noted  earlier,  the  primary  emphasis  of  this  part 
of  the  ISEG  evaluation  study  was  to  establish  the  compu¬ 
tational  capabilities/shortcomings  of  the  numerical  algori¬ 
thms  inherent  to  ADINA  (1977).  For  this  study,  it  was 
achieved  by  utilizing  specialized  benchmark  problems  with 
the  two-fold  purpose  of; 

i)  Performing  code  check-out;  and  more  importantly 

ii)  Establishing  the  algorithmic  sensitivities,  con¬ 
vergence  characteristics  and  typical  solution 
failure  modes  generic  to  the  ADINA  algorithmic 
apparatus. 

Since  ADINA  has  the  capability  to  handle  both  kinematic  and 
material  nonlinearity,  special  emphasis  was  given  to  ascer¬ 
taining  the  operating  characteristics  in  the  presence  of 
softening,  hardening,  elastic/plastic  and  load/unload  situa¬ 
tions. 

In  this  context,  the  benchmarking  served  to  establish 
the  sens: tivity/pathological  behavior  of  the  various  non¬ 
linear  algorithms  to  changes  in: 

i)  Convergence  criteria; 
ii)  Time  step  size; 

iii)  Load  increment  size; 

iv)  Material  models; 

v)  Order  of  integration  of  element  stiffness,  etc. 

To  enable  the  evaluation  of  the  various  behavioral  patho¬ 
logies,  specialized  coding  was  programmed  into  ADINA  to  allow 
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the  strain  energy,  potential  energy  and  kinetic  energy 
histories  to  be  evaluated  during  the  various  stages  of 
calculation.  This  allowed  for  the  monitoring  of  the  initia¬ 
tion  characteristics  of  anomalous  pathological  solution  be¬ 
havior. 

The  typical  types  of  solution  stops  encountered  during 
the  degradation  process  most  often  involved  either: 

i)  Out-of-balance  loads; 

ii)  Negative  pivots  or; 

iii)  Iteration  limit  stops; 

All  other  types  of  stops  (zero  Jacobian  etc.)  could  be  attri¬ 
buted  most  often  to  input  errors  in  the  data  preparation  or  to 
the  use  of  the  code  outside  its  range  of  capabilities.  For 
instance,  application  of  the  code  to  bifurcation  problems  in¬ 
volving  changes  from  positive  to  negative  definite  tangent 
stiffness  are  outside  the  range  of  application  of  the  solu¬ 
tion  algorithm.  Hence,  care  must  be  taken  in  dealing  with 
problems  involving  statically  or  dynamically  generated  post- 
buckling  behavior. 

Typically  it  was  found  that  the  nature  of  degradation 
is  dependent  on  the  relative  magnitude  of  load/time  incre¬ 
ment  size.  Usually  three  basic  types  of  solution  pathologies 
were  encountered.  Specifically: 

i)  Immediate  and  strong  nonmonotonicity; 

ii)  Moderate  but  progressively  increasing  non¬ 
monotonicity  and  nonpositive  definiteness;  and 

iii)  Mild  monotonicity  with  either  very  gradual  increases 
or  decreases  in  solution  oscillation. 


The  first  two  types  were  found  to  give  rise  to  out-of -balance 
loads  and  negative  pivots  (bifurcation  problems)  while  the 
third  type  caused  iteration  limit  stops. 

As  noted  earlier,  such  anomalous  behavior  is  not  an 
outgrowth  of  any  shortcomings  in  the  ADINA  architecture  and 
coding  practice.  Rather,  they  are  due  to  the  intrinsic  pro¬ 
perties  of  the  modified  NR  procedure  used  to  solve  the  non¬ 
linear  -field  equations  arising  in  both  the  static  and 
dynamic  branches.  Furthermore,  as  an  outgrowth  of  fairly  ex¬ 
tensive  code  check-out  benchmarking,  apart  from  minor  bugs, 
the  algorithmic  apparatus  inherent  to  ADINA  (1977)  is  both 
honestly  and  accurately  programmed. 

Due  to  the  results  of  the  sensitivity  studies,  it  has 
been  observed  that  while  ADINA  can  generally  be  forced  to 
generate  a  solution  to  highly  nonlinear  problems  via  the  vari¬ 
ous  restart  and  time  block  options,  such  an  approach  is  typi¬ 
cally  expensive  and  time  consuming.  Because  of  this  and  the 
adaptability  of  the  coding  and  architecture,  some  heuristic 
type  algorithmic  capabilities  should  be  incorporated  in  ADINA 
to  streamline  the  process  of  solution  generation  in  highly 
nonlinear  problems. 
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APPENDIX  1 :  Energy  Calculations 

As  was  first  noted  in  Chapter  III,  the  local  and  global 
strain  energy  are  used  to  monitor  the  computational  capa¬ 
bilities  and  shortcomings  of  the  algorithms  inherent  to  ADINA 
(1977).  For  the  present  purposes,  the  TL  version  of  the  energy 
functionals  will  be  derived.  The  development  will  include  both 
the  local  element  and  global  versions. 

To  start,  we  turn  to  the  NR  algorithm  defined  by 

[K_  ( Y) JAY  =  AR  (A-l) 

1  'V  %  'V. 

where  here  AY  denotes  the  increment  of  displacement  associated 
with  a  given  unbalanced  load  AR  and  [K„]  is  the  current  tangent 

%  A 

stiffness  matrix.  For  a  given  load  step,  the  iterative  solu¬ 
tion  of  (A-l)  obviously  leads  to  the  overall  incremental  res¬ 
ponse.  Graphically,  the  results  of  the  iterative  process  is 
defined  in  the  scalar  analogy  given  in  Fig.  Al-1.  The  incre¬ 
mental  energy  stored  during  a  given  iteration  step  is  es¬ 
sentially  the  shaded  area  illustrated  in  this  figure.  Reali¬ 
zing  that  the  ordinate  values  of  the  true  solution  surface  are 
given  by  the  relation ^ 

F  =  /  [B  (Y)]Ta(Y)dV  { A-2 ) 

^  J  %  %  % 

V 

it  follows  that  the  incremental  energy  stored  during  the 
iteration  step  can  be  approximated  by  the  inner  product 
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where  (  )  denotes  matrix  transposition,  and  B  is  the  element 
strain-nodal  displacement  transformation  matrix. 

In  terms  of  Fig.  (Al-2) ,  assuming  that  a  total  of  K1 
iteration  steps  are  associated  with  the  it^1  load  step,  then 
the  following  expression  can  be  developed  for  the  energy 
stored  namely 


K1 

I 

k=  1 


Ei 


-i  i 


k=  1 


Fk+1)T 


( A  -  4  ) 


Summing  (A-4)  over  the  entire  set  of  I  load  steps  yields  the 
requisite  overall  strain  energy  stored  for  a  given  problem 
namely 
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To  obtain  the  element  strain  energies  from  (A-5) ,  the 
various  nodal  displacement  and  force  vectors  must  be  inter¬ 
preted  from  a  local  point  of  view  namely 


k  k+1  k 

where  here  Y.  Y.  and  AY.  denote  local  element  partitions 

Opie  'vie 

Jr  k+1  k 

of  YV,  Y.  and  AY.  and  R  is  the  region  occupied  by  the 
<vi  %i  %i  o  e  r  r  i 

element. 

Because  of  the  form  of  either  (A-5)  or  (A-6) ,  any  form 
of  tangent  stiffness  type  constitutive  law  can  be  accomodated. 
This  is  particularly  true  for  those  inherent  to  ADINA  (1977). 


.  _ , 
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The  actual  calculations  programmed  into  the  ADINA' 
algorithmic  flow  are  depicted  in  Appendix  3.  As  noted 
earlier,  such  options  include  the  capability  of  handling 

1)  Multiple  element  groups; 

2)  In/out  of  core  solution  options; 

3)  All  available  material  property  options  and; 

4)  All  algorithmic  options  (static  and  dynamic) . 
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APPENDIX  2  :  Code  Checkout 

The  various  benchmarks  that  were  run  strictly  for 
preliminary  code  check  out  will  be  briefly  discussed  in 
this  Appendix.  The  problem  types  employed  fell  into  three 
main  structural  categories  namely: 

1)  Beam/ truss  type  structure; 

2)  2-D  continuum  and; 

3)  3-D  continuum. 

The  analyses  ran  consisted  of  static,  dynamic  and  eigenvalue/ 
vector  problems. 

For  check  out  pruposes,  NFAP ,  NONSAP,  STRUDL  and  NASTRAN* 
were  employed  to  perform  parallel  runs.  This  enabled  both 
linear  and  nonlinear  verifications  to  be  performed  on  well  tried 
and  readily  available  "general"  purpose  codes. 

Note  the  linear  benchmarks  with  STRUDL  and  NASTRAN  en¬ 
abled  a  thorough  determination  of  any  anomalies  in  much  of  the 
ADINA  bookkeeping  associated  with: 

1.  In/out  of  core  solutions; 

2.  Handling  of  various  boundary  conditions; 

3.  Handling  of  different; 

i)  element  groups 

ii)  material  properties 

iii)  node  sequencing  etc. 

Additionally,  such  tests  permitted  the  benchmarking  of  the 
initial  structure  stiffness. 


*Cosmic  version  levels  16.  and  17.5 
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Once  such  factors  were  verified,  the  various  geometric 
and  material  nonlinearities  were  benchmarked  within  the  capa¬ 
bilities  of  NONSAP  and  its  derivative  NFAP . 

The  primary  source  of  benchmarks  was  drawn  from  the 
ADINA  sample  problems.  Whenever  possible,  their  results  were 
verified  via  either  NONSAP,  NFAP  or  the  linear  initial  stiff¬ 
ness  results  via  STRUDL  and  NASTRAN.  The  benchmarks  con¬ 
sisted  of  the  following  list  of  problems  namely: 

1.  Tower  Cable; 

i)  static 

ii)  frequency 

2.  Rubber  sheet;  Mooney  Rivlin  material,  small/ 
large  deformation 

i)  static 

ii)  dynamic 

iii)  frequency  (linear) 

3.  Spherical  shell;  elastic,  plastic,  concrete, 
small/large  deformation 

i)  static 
ii)  dynamic 

iii)  frequency  (linear) 

iv)  onset  of  static/dynamic  buckling 

4.  Thick  walled  cylinder;  thermo-elastic-plastic, 
small/large  deformation 

i)  static 
ii)  dynamic 


5.  Simply  supported  plate;  elastic,  plastic,  concrete 
small/large  deformation 

i)  static 

ii)  dynamic 

iii)  frequency/linear 

6.  Arch;  elastic ,  plastic,  concrete,  small/large 
deformation 

i)  static 

ii)  dynamic 

iii)  frequency  (linear) 

iv)  onset  of  static/dynamic  buckling 

7.  Cyclic  creep  analysis  of  a  thick  walled  cylinder. 

8.  Underground  opening,  elastic,  plastic,  concrete, 
static  analysis. 

9 .  Reinforced  concrete  beam 

i)  static 

ii)  dynamic 

iii)  frequency 

10.  Pipe  whip;  dynamic  response. 

During  the  initial  attempts  to  run  the  foregoing  pro¬ 
blems,  benchmarks  1-7  yielded  successful  solutions  for  all 
the  various  analyses  performed  (static,  dynamic,  frequency, 
etc.).  Such  was  not  the  case  for  problems  8-10 .  In  par¬ 
ticular,  in  our  version  of  ADINA  (1977)  ,  negative  pivots 
were  initially  encountered.  After  some  debugging,  dynamic 
allocation  problems  were  found  in  the  MAXA  array.  Once 


corrected,  all  the  above  listed  problems  yielded  successful 
solutions  for  the  various  analyses  performed. 

To  check  the  various  bookkeeping  features  of  ADINA 
alluded  to  earlier,  (in/out  of  core  solutions,  etc.),  the 
foregoing  problems  were  modified  in  several  different  ways 
namely; 

1.  Elementing  increased/decreased 

2.  Element  type  changed 

3.  Node  sequencing  changed 

4.  Changes  introduced  in 

i)  Material  types 

ii)  Boundary  conditions 

iii)  Loadings  etc. 

In  order  to  verify  the  in  and  out  of  core  features,  the 
dynamic  allocation  and  block  size  was  increased  and  decreased. 
Such  an  approach  allowed  the  same  problem  to  be  run  in  and 
out  of  core.  For  all  the  foregoing  manipulations,  apart  from 
the  minor  allocation  difficulties  associated  with  the  MAXA 
array,  no  difficulties  or  aberration  were  encountered  in  either 
the  beam,  2-D  or  3-D  elements. 

In  addition  to  the  foregoing  problems,  several  additional 
code  check  out  benchmarks  were  considered.  The  main  thrust 
of  this  group  of  problems  was  to  establish  the  capability  of 
ADINA  to  predict  the  onset  of  static  and  dynamic  structural 
instability.  In  this  context,  buckling  problems  involving 
beams,  arches,  rings  and  spherical  caps  were  considered.  This 


A2-5 


included: 

1.  Offset  loaded  beam/columns  with  varying  boundary 
conditions ; 

2.  Arches  with  varying  depths  and  boundary  conditions; 

3.  Radially  loaded  rings  and; 

4.  Spherical  caps  with  varying  depths  and  boundary 
conditions . 

The  results  of  the  above  noted  "stability"  problems 
were  compared  with  NONSAP  and  NFAP  for  independent  veri¬ 
fication.  Note  within  the  limits  of  the  tangent  modulus/MNR 
approach,  good  resolution  was  obtained  for  such  problems. 

Such  was  not  the  case  though  for  the  beam  element.  In  par¬ 
ticular,  for  offset  column  type  buckling  problem,  immediate 
solution  failure  was  encountered.  From  the  nature  of  the  pro¬ 
gram  stop,  (negative  pivot) ,  the  difficulty  appears  to  be 
isolated  to  the  beam  element.  This  also  follows  from  the  fact 
that  good  buckling  predictions  were  obtained  via  the  2-D  TL 
and  UL  formulation  of  the  same  problem.  Apart  from  the  fore¬ 
going  difficulties,  ADINA  (1977)  appears  to  be  largely  error 


free. 


1)  Static  Analysis  of  a  Rubber  Sheet 

2)  Dynamic  Analysis  of  a  Spherical  Cap 

3)  Dynamically  Loaded  Cantilever  Beam 
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